#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <stdlib.h>
#define ll long long
#define maxn 100+5
#pragma GCC optimize(2)
const double eps=1e-11;

int QR(double a[][maxn][maxn],double Q[][maxn][maxn],double R[][maxn],double c[],int n);
double sgn(double x);
int allZero(double a[][maxn],int r,int n);
void rev(double a[][maxn],double b[][maxn],int n);

int main(){
    double a[2][maxn][maxn],Q[2][maxn][maxn],R[maxn][maxn],c[maxn];
    int n;
    printf("输入矩阵的阶数n:");
    scanf("%d",&n);
    printf("输入%d*%d矩阵:\n",n,n);
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
            scanf("%lf",&a[0][i][j]);
    memset(Q,0,sizeof(Q));
    memset(R,0,sizeof(R));
    int now=QR(a,Q,R,c,n);
    printf("Q=\n");
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n;++j)
            printf("%10.6lf",Q[now][i][j]);
        printf("\n");
    }
    printf("\nR=\n");
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n;++j){
            printf("%10.6lf",a[now][i][j]);
        }
        printf("\n");
    }
    return 0;
}

int QR(double a[][maxn][maxn],double Q[][maxn][maxn],double R[][maxn],double c[],int n){
    int now=1,pre=0;
    //double c[maxn];
    for(int i=1;i<=n;++i) Q[0][i][i]=1;
    for(int r=1;r<=n-1;++r){
        if(allZero(a[pre],r,n)){
            for(int i=1;i<=n;++i)
                for(int j=1;j<=n;++j)
                    Q[now][i][j]=Q[pre][i][j],a[now][i][j]=a[pre][i][j];
            continue;
        }
        double d=0.0,h;
        for(int i=r;i<=n;++i)
            d+=(a[pre][i][r]*a[pre][i][r]);
        d=sqrt(d);
        c[r]=-1.0*sgn(a[pre][r][r])*d;
        h=c[r]*c[r]-c[r]*a[pre][r][r];
        double u[maxn];
        memset(u,0,sizeof(u));
        u[r]=a[pre][r][r]-c[r];
        for(int i=r+1;i<=n;++i){
            u[i]=a[pre][i][r];
        }

        double omega[maxn];// omega=Q[r]*u[r]
        for(int i=1;i<=n;++i){
            double sum=0.0;
            for(int j=1;j<=n;++j)
                sum+=Q[pre][i][j]*u[j];
            omega[i]=sum;
        }
        for(int i=1;i<=n;++i){//Q[r+1]=Q[r]-omega*u^t /h
            for(int j=1;j<=n;++j)
                Q[now][i][j]=Q[pre][i][j]-omega[i]*u[j]/(h*1.00);
        }

        double reverseA[maxn][maxn],p[maxn];
        rev(a[pre],reverseA,n);
        for(int i=1;i<=n;++i){
            double sum=0.0;
            for(int j=1;j<=n;++j)
                sum+=reverseA[i][j]*u[j];
            p[i]=sum/h;
        }

        for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j)
                a[now][i][j]=a[pre][i][j]-u[i]*p[j];

        int temp=now;
        now=pre;
        pre=temp;
    }
    return now==1?0:1;
}

inline void rev(double a[][maxn],double b[][maxn],int n){
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
            b[i][j]=a[j][i];
}

int allZero(double a[][maxn],int r,int n){
    for(int i=r+1;i<=n;++i)
        if(a[i][r]!=0)
            return 0;
    return 1;
}

inline double sgn(double x){
    return (x-0.0<eps?-1.0:1.0);
}